// System includes
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <time.h>
using namespace std;

// Matrix configurations
#define N 4                                           // Matrix order
#define BLOCK_SIZE 2                                  // Block size; N % BLOCK_SIZE must = 0

#define N_SQUARED (N * N)                             // Total number of cells
#define MAX_RAND 10                                   // Random range upper bound for populating A and B
#define BLOCK_SIZE_SQUARED (BLOCK_SIZE * BLOCK_SIZE)  // Pre-calculated to save a few multiplications

// Matrices
int* A;
int* B;
int* C;

// Forward function declarations
void multiply();
void multiply_block();
void block_add(int i_block, int j_block, int k_block);
void print_matrix(const char* name, int* m);

//
// Utility function to print a matrix to standard output.
//
void print_matrix(const char* name, int* m)
{
   std::cout << name << "=[" << std::endl;
   for (int i = 0; i < N; ++i)
   {
      std::cout << "\t";
      for (int j = 0; j < N; ++j)
      {
         std::cout << m[i * N + j] << "\t";
      }
      std::cout << std::endl;
   }
   std::cout << "]" << std::endl;
}

//
// Multiply two matrices together using standard matrix multiplication.
//
void multiply() 
{
   for (int i = 0; i < N; ++i) 
   {
      for (int j = 0; j < N; ++j) 
      {
         C[i * N + j] = 0;
         for (int k = 0; k < N; ++k)
         {
            C[i * N + j] += A[i * N + k] * B[k * N + j];
         }
      }
   }
}

//
// Perform matrix multiplication on only a part of a larger matrix.
// The block being calculated is indexed by the provided parameters.
//
void block_add(int i_block, int j_block, int k_block)
{
   // Load sub/block matrices from A, B, C of size BLOCK_SIZExBLOCK_SIZE
   int* a = new int[BLOCK_SIZE_SQUARED];
   int* b = new int[BLOCK_SIZE_SQUARED];
   int* c = new int[BLOCK_SIZE_SQUARED];
   for (int i = 0; i < BLOCK_SIZE; ++i)
   {
      for (int j = 0; j < BLOCK_SIZE; ++j)
      {
         for (int k = 0; k < BLOCK_SIZE; ++k)
         {
            a[i * BLOCK_SIZE + k] = A[(i_block + i) * N + (k_block + k)];
            b[k * BLOCK_SIZE + j] = B[(k_block + k) * N + (j_block + j)];
            c[i * BLOCK_SIZE + j] = C[(i_block + i) * N + (j_block + j)];
         }
      }
   }

   // Multiply just the blocks
   for (int i = 0; i < BLOCK_SIZE; ++i)
      for (int j = 0; j < BLOCK_SIZE; ++j)
         for (int k = 0; k < BLOCK_SIZE; ++k)
            c[i * BLOCK_SIZE + j] += (a[i * BLOCK_SIZE + k] * b[k * BLOCK_SIZE + j]);

   // Store block that we computed back into the whole matrix
   for (int i = 0; i < BLOCK_SIZE; ++i)
      for (int j = 0; j < BLOCK_SIZE; ++j)
         for (int k = 0; k < BLOCK_SIZE; ++k)
            C[(i_block + i) * N + (j_block + j)] = c[i * BLOCK_SIZE + j];

   // Free memory
   delete [] a;
   delete [] b;
   delete [] c;
}

//
// Multiply two matrices together using block matrix multiplication.
// This method divides the entire matrix up into equal sized sub-matrices
// called blocks and adds up the results in pieces.  
//
void multiply_block()
{
   // Zero out C; this can be done prior to this function if need be.
   // I only do it here because I run the normal multiply() for comparison
   // purposes.
   for (int i = 0; i < N_SQUARED; ++i)
      C[i] = 0;

   // Multiply sub matrices
   for (int i_block = 0; i_block < N; i_block+=BLOCK_SIZE)
      for (int j_block = 0; j_block < N; j_block+=BLOCK_SIZE)
         for (int k_block = 0; k_block < N; k_block+=BLOCK_SIZE)
            block_add(i_block, j_block, k_block);
}

//
// Main entry point for program.
//
int main()
{
   // Initialize random seed
   srand(time(NULL));

   // Allocate memory for matrices
   A = new int[N_SQUARED];
   B = new int[N_SQUARED];
   C = new int[N_SQUARED];

   // Initialize A and B random elements
   for (int i = 0; i < N; ++i)
   {
      for (int j = 0; j < N; ++j)
      {
         A[i * N + j] = rand() % MAX_RAND + 1;
         B[i * N + j] = rand() % MAX_RAND + 1;
      }
   }
   
   // Print matrices A and B
   //print_matrix("A", A);
   //print_matrix("B", B);
   
   // Multiply the matrices normally
   multiply();
   print_matrix("C=AxB", C);

   // Multiply the matrices using block decomposition
   multiply_block();
   print_matrix("C=AxB", C);

   // Free memory
   delete [] A;
   delete [] B;
   delete [] C;
}